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ABSTRACT 

The Fermi Large Area Telescope (LAT) reported the first definitive GeV detections of the binaries LS 1 +61°303 
and LS 5039 in the first year after its launch in 2008 June. These detections were unambiguous as a consequence 
of the reduced positional uncertainty and the detection of modulated y-ray emission on the corresponding orbital 
periods. An analysis of new data from the LAT, comprising 30 months of observations, identifies a change in the 
y-ray behavior of LS I +61° 303. An increase in flux is detected in 2009 March and a steady decline in the orbital 
flux modulation is observed. Significant emission up to 30 GeV is detected by the LAT; prior data sets led to upper 
limits only. Contemporaneous TeV observations no longer detected the source, or found it — in one orbit — close to 
periastron, far from the phases at which the source previously appeared at TeV energies. The detailed numerical 
simulations and models that exist within the literature do not predict or explain many of these features now observed 
at GeV and TeV energies. New ideas and models are needed to fully explain and understand this behavior. A detailed 
phase-resolved analysis of the spectral characterization of LS I +61° 303 in the GeV regime ascribes a power law 
with an exponential cutoff spectrum along each analyzed portion of the system’s orbit. The on-source exposure of 
LS 5039 is also substantially increased with respect to our prior publication. In this case, whereas the general y-ray 
properties remain consistent, the increased statistics of the current data set allows for a deeper investigation of its 
orbital and spectral evolution. 
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(LS I +61 303) - X-rays: individual (LS 5039) 
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1. INTRODUCTION 

To date there are only a handful of X-ray binaries that have 
been detected at high (HE; 0. 1-100 GeV) or very high energies 

(>100 GeV): LS I +61°303 (Albert et al. 2006; Acciari et al. 

2008; Abdo et al. 2009a), LS 5039 (Aharonian et al. 2005b; 

Abdo et al. 2009b), PSR B 1259— 63 (Aharonian et al. 2005a; 

Abdo et al. 2011; Tam et al. 2011), Cyg X-3 (Abdo et al. 2009c), 

and Cyg X— 1 (Albert et al. 2007; Sabatini et al. 2010). Recently, 

two new binaries were found: 1FGL J1018.6 — 5856, with a pe- 
riod of 16.6 days found in the GeV regime (Corbet et al. 2012) 
and HESS J0632+057 (Falcone et al. 2011; Ong 2011; 

Mariotti 2011), for which a period of ~320 days was detected in 
X-rays (Bongiorno et al. 2011). Of these sources only 
LS I +61°303, LS 5039, and PSR B 1259-63 share the property 
of being binaries detected at both GeV and TeV energies. The 
other systems have been unambiguously detected only in one 
band, either at GeV or at TeV, see, e.g., the case of Cyg X— 3 in 
Aleksic et al. (2010). In the case of Cyg X— 1, with the hint of 
TeV detection itself being at the level of four standard deviations 

(4a), claims of detection at GeV energies by the Astrorivelatore 

Gamma a Immagini Leggero remain uncertain with concurrent 


Fermi Large Area Telescope (LAT) observations (Hill et al. 
2011). It is yet uncertain whether these spectral energy distribu- 
tion (SED) differences reflect an underlying distinct nature, or 
are just a variability signature in different bands. 

The nature of the binary compact object in LS I +61°303, 
LS 5039, HESS J0632+057, and 1FGL J1018.6— 5856 is as 
yet undetermined (Hill et al. 2011). Both neutron star (e.g., 
PSR B 1259— 63) and probable black hole (e.g., Cyg X— 3) 
binary systems have been detected at GeV energies, and so 
both types of compact object are viable in the undetermined 
systems. Recently, the Burst Alert Telescope on board Swift 
reported a magnetar-like event which may have emanated from 
LS I +61°303 (Barthelmy et al. 2008; Torres et al. 2012). If true, 
then this would be the first magnetar found in a binary system. 

The early LAT reports of GeV emission from LS 5039 and 
LS I +61° 303 were based upon 6-9 months of survey observa- 
tions (Abdo et al. 2009a, 2009b). Both sources were detected at 
high significance and were unambiguously identified with the 
binaries by their flux modulation at the corresponding orbital 
periods, 26.4960 days for LS I +61°303 (Gregory 2002) and 
3.90603 days for LS 5039 (Casares et al. 2005). The modula- 
tion patterns were roughly consistent with expectations from 
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inverse Compton scattering plus y-y absorption models, and 
were anti-correlated in phase with pre-existing TeV measure- 
ments (e.g., Albert et al. 2009; Aharonian et al. 2006). The 
anti-correlation of GeV-TeV fluxes is, in fact, a generic feature 
of these models, where the GeV emission is enhanced (reduced) 
when the highly relativistic electrons seen by the observer en- 
counter the seed photons head-on (rear-on) (e.g., see Bednarek 
2007; Sierpowska-Bartosik & Torres 2007; Dubus et al. 2008; 
Khangulyan et al. 2008). Fermi-L AT measurements provided a 
confirmation of these predictions. 

The spectra of both sources were best modeled with exponen- 
tial cutoffs in their HE spectra, at least along part of the orbit. 
Specifically, an exponential cutoff was statistically a better fit to 
the SED compared with a pure power law at phases surround- 
ing the superior conjunction (SUPC) of LS 5039 and in the 
orbitally averaged spectrum of LS I +61°303. Statistical limita- 
tions of the data prevented the determination or the ruling out of 
an exponential cutoff in any part of the orbit of LS I +61° 303 or 
in the inferior conjunction (INFC) of LS 5039. The SEDs with 
the exponential cutoffs that were reported were reminiscent of 
the many pulsars the LAT has detected since launch (Abdo et al. 
2009b), although this was far from a proof of their pulsar nature. 
To date no pulsations have been found at GeV energies, or at any 
other wavelengths, despite deep dedicated searches (see, e.g., 
Rea et al. 2010, 2011). 

Since Fermi was launched, both the Major Atmospheric 
Gamma-ray Imaging Cherenkov Telescopes (MAGIC) and the 
Very Energetic Radiation Imaging Telescope Array (VERITAS) 
have performed observations of LS I +61°303. No TeV detection 
was reported after 2008 October, until the source unexpectedly 
reappeared, once, at periastron (Acciari et al. 2011; Ong 2010). 
At the same time, a hard X-ray multi-year analysis (Zhang et al. 
2010) and a long-term X-ray campaign on LS I +61°303 using 
the Rossi X-ray Timing Explorer (RXTE) has been conducted 
covering the whole extent of the LAT observations (see Torres 
et al. 2010; Li et al. 2011). In addition, simultaneous and archival 
data from long-term monitoring of radio and Ha emission is 
available for comparison in a multi-wavelength context. In this 
work, we present the results of the analysis of 30 months of LAT 
survey observations of both LS I +61°303 and LS 5039. We 
investigate the long-term flux variations of the sources, as well 
as variations in the amplitude of their orbital flux modulation, 
and we explore the possible spectral variability for both systems, 
finally putting and interpreting these observations in the context 
of the source behavior at other frequencies. 

2. OBSERVATIONS AND DATA REDUCTION 

The Fermi- LAT is an electron-positron pair production tele- 
scope, featuring solid state silicon trackers and cesium iodide 
calorimeters, sensitive to photons from ~20 MeV to >300 GeV 
(Atwood et al. 2009). It has a large ~2.4 sr field of view (at 
1 GeV) and an effective area of ~8000 cm 2 for > 1 GeV. 


2.1. Data Set 

The Fermi survey mode operations began on 2008 August 4; 
in this mode, the observatory is rocked north and south on 
alternate orbits to provide a more uniform coverage, so that 
every part of the sky is usually observed for ~30 minutes every 
~3 hr. Therefore, the two sources of interest were monitored 
regularly without significant breaks, allowing us to draw a 
complete picture of their behavior in y-rays over the last 


2 years. The data set used for this analysis spans 2008 August 4 
through 201 1 January 24. 

The data were reduced and analyzed using the Fermi Science 
Tools v9r20 package. 12 The standard on board filtering, 
event reconstruction, and classification were applied to the 
data (Atwood et al. 2009). The high-quality “diffuse” event 
class was used together with the “Pass 6 v3 Diffuse” instru- 
ment response functions (IRFs). Time periods when the tar- 
get source was observed at a zenith angle greater than 105° 
were excluded to limit contamination from Earth limb pho- 
tons. Where required in the analysis, models for the Galactic 
diffuse emission ( gll_iem_v02.fit ) and isotropic backgrounds 
(isotropic_iem_v02.txt) were used. 13 


2.2. Spectra! Analysis Methods 

The binned maximum likelihood method of gtlike, included 
in the Science Tools, was used to determine the intensities and 
spectral parameters presented in this paper. We used all photons 
with energy >100 MeV in a circular region of interest (ROI) 
of 10° radius centered at the position of LS I +61°303 and 
LS 5039, respectively. For source modeling, the 1FGL catalog 
(Abdo et al. 2010a), derived from 1 1 months of survey data, and 
the first Fermi pulsar catalog (Abdo et al. 2010b) were used; all 
sources within 15° of the ROI center were included. The energy 
spectra of point sources included in the catalog within our ROI 
are modeled by a simple power law. 


dN 

d~E 



(1) 


with the exception of known y-ray pulsars, which were modeled 
by power laws with exponential cutoffs described by 


dN 

Je 



E 

Teuton 


(2) 


The spectral parameters were fixed to the catalog values, except 
for the sources within 3 deg of the candidate location. For 
these latter sources, the flux normalization was left free. All 
of the spectral parameters of the two subject binaries were left 
free for the fit. Source detection significance is determined using 
the Test Statistic value, TS = —2 ln(Lo/L i ) which compares 
the likelihood ratio of models including, e.g., an additional 
source, with the null hypothesis of background only (Mattox 
et al. 1996). 

To estimate the systematic errors, which are mainly caused 
by uncertainties in the effective area and energy response of 
the LAT as well as background contamination, we use the so- 
called bracketing IRFs. These are IRFs with effective areas that 
bracket those of our nominal IRF above and below by linearly 
connecting differences of (10%, 5%, 20%) at log(E/MeV) of 
(2, 2.75, 4), respectively. 


2.3. Timing Analysis Methods 

Light curves are extracted using aperture photometry, taking 
an aperture radius of 1° and using the gtbin tool. The exposure 
correction is performed with the tool gtexposure assuming 
the spectral shape of the source to be a power law with 


12 See the Fermi Space Science Center (FSSC) Web site for details of the 
Science Tools: http://fermi.gsfc.nasa.gov/ssc/data/analysis/. 

13 Descriptions of the models are available from the FSSC: 
http://fermi.gsfc.nasa.gov/ssc. 
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Figure 1 . Left: spectral energy distribution of LS 5039 as observed by the Fermi-LAT (red). The gray open and filled circles are H.E.S.S. spectra in INFC and SUPC, 
respectively (Aharonian et al. 2006) and are not simultaneous with the GeV measurements. Right: Fermi - LAT and H.E.S.S. spectra of LS 5039 in SUPC (blue) and 
INFC (red). 

(A color version of this figure is available in the online journal.) 


an exponential cutoff (see Sections 3.1 and 4.1). These light 
curves are not background subtracted. The folded light curves 
shown in the subsequent sections are derived by performing 
gtlike fits for each phase bin. Therefore, all of them are 
effectively background subtracted. We check that both methods 
for generating light curves, aperture photometry, and gtlike 
fits are consistent with each other when the former light curves 
are background subtracted too. 

The primary method of timing analysis employed searches for 
periodic modulation by calculating the weighted periodogram 
of the light curve (Lomb 1976; Scargle 1982; Corbet & Dubois 
2007). The light curve is constructed by summing, for each 
photon, the estimated probability that the photon came from 
the source of interest. The probability will be both spatially 
and spectrally dependent. Because this technique allows for the 
correct weighting of each photon, it intrinsically improves the 
signal-to-noise and allows the use of a larger aperture. This 
method has successfully been applied to increase the LAT 
sensitivity for the detection of pulsars (Kerr 2011). However, 
in the basic form of this technique, the weight for any particular 
energy /position is fixed. This means that changes in source 
brightness will not be reflected in the weights and can result 
in incorrect probabilities. The calculation of probabilities was 
performed using the tool gtsrcprob and the same source 
model file derived from the 1FGL catalog and used in the 
spectral analysis. Since the exposure of the time bins was 
variable, the contribution of each time bin to the power spectrum 
was weighted based on its relative exposure. Period errors are 
calculated using the method of Horne & Baliunas (1986). 

3. LS 5039 RESULTS 

LS 5039 is located in a complicated region toward the 
inner Galaxy with high Galactic diffuse emission and many 
surrounding y-ray sources. In particular, the LAT detected a 
bright (8.7 x 10~ 7 photons cm -2 s -1 above 100 MeV) y-ray 
pulsar, PSR J1826-1256, ~2° away from LS 5039. Following 
the analysis performed in the earlier LAT paper (Abdo et al. 


2009b), we discarded events whose arrival times correspond 
to the peaks of the pulsar cycle of PSR J 1826— 1256 in order 
to minimize the contamination from the pulsar. The excluded 
pulse phase of PSR J1826— 1256 is 0.05 < <p p < 0.2 and 
0.625 < (fi p < 0.775 (see Figure 37 in Ray et al. 2011), which 
results in a loss of 30% exposure on LS 5039. To account for the 
loss, a scaling factor of 1 /0.7 is multiplied to fluxes obtained 
with maximum likelihood fits. 

3.1. Orbitally Averaged Spectrum 

The orbitally averaged spectrum of LS 5039 was initially 
investigated by fitting a power law and a power law with an 
exponential cutoff to the data. We compare two models utilizing 
the likelihood ratio test (Mattox et al. 1996), i.e., for the ratio 
2 x Alog(Likelihood), we assume a y 2 -distribution to calculate 
the probabilities taking into account the corresponding degrees 
of freedom (Eadie et al. 1971). In this case, the significance 
of a spectral cutoff was assessed by comparing the likelihood 
ratio between the power law and cutoff power-law cases, which 
is —2 ln(LpL/L C utoff) = 94.9, where Lpl and L cut off are the 
likelihood values obtained for the spectral fits with a power 
law and a cutoff power law, respectively. This indicates that 
the simple power-law model is rejected at the 9.7er level in 
favor of the cutoff power law. The best-fit parameters for the 
cutoff power-law model are T = 2.06 ± 0.06 sta t ± 0.11 syst 
and Ecutoff = 2.2 ± 0.3 stat ± 0.5 syst GeV with a flux of 
E> ioo MeV = (6.1 ± 0.3 stat ± 2.1 syst ) x 10~ 7 cm~ 2 s -1 integrated 
above 100 MeV. Using a cutoff power-law spectral model, the 
maximum likelihood fit yields a test statistic of TS = 1623 for 
the LS 5039 detection; equivalent to ~40cr. We also tried a 
broken power-law spectral model for LS 5039 in addition to 
an exponentially cutoff power law. We found that the broken 
power law gives lower T S values than the exponentially cutoff 
power-law case. 

Spectral points in each energy band were obtained by dividing 
the data set into separate energy bins and performing maximum 
likelihood fits for each of them. The resulting SED is plotted 
in Figure 1 together with the best-fit cutoff power-law model. 
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Figure 2. Light curve of LS 5039 in time bins of six orbital cycles between 
100 MeV and 300 GeV. One orbit is 3.90532 day long. 

(A color version of this figure is available in the online journal.) 


Interestingly, the SED shows significantly higher flux (one 
spectral data point) at >10 GeV than the expected flux from 
the best-fit cutoff power law, possibly suggesting another 
component at HEs. 

One idea explored for these sources (especially, for LS 5039, 
see, e.g., Torres 2011) is that the y -ray emission could be 
understood as having two components: one would be the 
magnetospheric GeV emission from a putative pulsar and the 
other from the inter- wind region or from the pulsar wind zone. 
The latter would be unpulsed and would vary with the orbital 
phase, the former would be steady and pulsed. The current data 
for LS 5039 would indeed allow for this possibility, especially 
because of the possible HE component found in the GeV 
spectrum. 

To test the significance of the additional component, we added 
to the model a power-law source at the location of LS 5039 in 
addition to the cutoff power-law source. Model A is a power 
law with a cutoff and model B a power law with a cutoff plus 
an additional power law. Thus, model B has two more free 
parameters compared to model A. According to the likelihood 
ratio test, the probability of incorrectly rejecting model A is 
6.1 x 10 -6 (4.5er). 

The best-fit parameters for the putative additional component 
are T = 1.6 ± 0.4 stat ± 0.3 syst and E>i 0 GeV = (1-6 ± 0.6 stat ± 
0.7 syst ) x 1 0 10 cm 2 s' 1 . The addition of the HE component 
slightly affects the parameters of the cutoff power law. The best- 
fit parameters for the latter are T = 2.02 ± 0.06 stat ± 0.10 syst 
and Ecutoff = 2.0 ± 0.3 sta t ± 0.4 syst GeV with a flux of 
E>iooMeV = (6.0 ± 0.3 sta t ± 1.9 syst ) X 10” 7 cm' 2 s' 1 integrated 
above 100 MeV. 


3.2. Phase-resolved Analysis 

Following the H.E.S.S. analysis by Aharonian et al. (2006), 
as well as the previous one, the whole data set was divided into 
two orbital intervals: SUPC (<j) < 0.45 and 0.9 < <p) and INFC 
(0.45 < <p < 0.9). The SUPC and INFC data were analyzed in 
the same way as the orbitally averaged data. Being consistent 
with our previous paper, the power-law assumption for the SUPC 
spectrum is rejected with —2 ln(L PL /L cllt0 ff) = 81.2, or at a 
rejection significance of ~9er. The best-fit parameters are T = 
2.07 ± 0.07 stat ± 0.08 syst , Ecutoff = 1.9±0.3 sta t±0.3 sys t GeV, 
and ioo MeV = (7.8 ± 0.4 stat ± 1.9 syst ) x 10 _7 cm~ 2 s _1 . 

Although a single power law was not rejected for INFC in 
our previous analysis using 10 months of data (Abdo et al. 
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Figure 3. Lomb-Scargle power spectrum of the whole GeV data set on LS 5039. 
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2009b), a cutoff power law is preferred also for INFC with 
the present data set. The likelihood ratio for the INFC data 
is —2 ln(LpL/L C utoff) = 21.7, which corresponds to 4.7cr. 
The parameters for the INFC spectrum are T = 1.99 ± 
0.13 sta t ± 0.07 syst , Ecutoff = 2.6 ± 0.7 stat ± 0.9 syst GeV, and 
E>iooMeV = (3.9±0.4 stat ± 1.5 syst ) x 10' 7 cm' 2 s” 1 . Therefore, 
the SUPC and INFC spectral shapes are completely consistent 
with one another within the errors. The only difference is the 
normalization and hence the total flux. On the other hand, the 
spectrum for INFC (red points in the right panel of Figure 1) 
seems to exhibit additional structure below 1 GeV. The limited 
statistics and the large contribution of diffuse emission at low 
energies, however, prevents solid conclusions on whether a more 
complicated fit (e.g., a double broken power law or a broken 
power law with a cutoff) would be preferred. 

We also searched for emission from the HE component in 
the SUPC and INFC spectra. However, the TS of the additional 
components compared with a power law with exponential cutoff 
are only 13.6 and 10.9 for SUPC and INFC, respectively, 
and do not confirm a second spectral component. The SUPC 
and INFC SEDs were obtained using the same method as the 
orbitally averaged spectra and are plotted in the right panel of 
Figure 1. 


3.3. Light Curve 

Figure 2 shows the light curve for LS 5039 over 30 months de- 
rived by performing gtlike fits on time bins which contain six 
orbital cycles each. The light curve for LS 5039 does not show 
any significant flux changes. Constructing the periodogram of 
the weighted photon light curve yields a significant detection of 
a periodicity at 3.90532 ± 0.0008 days. This is consistent with 
the known orbital period of LS 5039. The Lomb-Scargle power 
spectrum of LS 5039 is shown in Figure 3. The stability of the 
orbital modulation was investigated and no significant variation 
in the modulation fraction as a function of time was found. 

4. LSI +61°303 RESULTS 
4.1. Orbitally Averaged Spectral Analysis 

We have derived the spectrum of orbitally averaged LAT 
data, i.e., without any selection criteria (cuts) concerning the 
orbital phase, for the LS I +61° 303 system. The spectral points 
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Table 1 

T S Values for LS I +61° 303 for Different Spectral Shapes (See Section 4.3 for Details) 


Data Set 

Power Law + Exponential Cutoff 

Power Law 

Broken Power Law 


TS 

TS 

TS 

30 months of data 

23995 

23475 

23970 

Data before 2009 March 

3404 

3314 

3415 

Data after 2009 March 

20714 

20283 

20699 

Inferior conjunction (geometrically) 

12548 

12326 

12512 

Superior conjunction (geometrically) 

11711 

11422 

11700 

Inferior conjunction (angle cut) 

6670 

6562 

6665 

Superior conjunction (angle cut) 

6083 

5986 

6063 

Periastron 

11656 

11450 

11636 

Apastron 

12377 

12059 

12361 


Table 2 



Parameters for LS I +61 c 

;> 303 from Spectral Fitting with Power Law with Exponential Cutoff (See Section 4.3 for Details) 

Data Set 

Photon Index V 

Cutoff Energy 

Flux >100 MeV 



(GeV) 

(x 10 -6 photons cm" 2 s _1 ) 

First 8 months of data 

2.21 ± 0.04 s tat ± 0.06 syst 

6.3 zt : 1 - 1 stat i 0-4-syst 

0.82 dz 0.03 sta t dz 0.07 S y S t 

30 months of data 

2.07 dz 0.02 sta t dz 0.09 sys t 

3.9 dz 0.2 s tat dz 0.7 S y S t 

0.95 ± 0.01 s tat ± 0.07 syst 

Data before 2009 March 

2.08 ± 0.04 s tat ± 0.09 sys t 

4.0 ± 0.6stat ± 0.7 syst 

0.75 ± 0.03 st at ± 0.07 syst 

Data after 2009 March 

2.07 ± 0.02 s tat ± 0.09 sys t 

3.9 dz 0.3 S tat dz 0.7 S yst 

1 .00 ± 0.01 s tat ± 0.07 syst 

Inferior conjunction (geometrically) 

2.14 ± 0.02 sta t ± 0.09 sys t 

4.0 dz 0.4 s tat dz 0.7 S y S t 

1 .07 ± 0.02 stat ± 0.07 S y St 

Superior conjunction (geometrically) 

2.02 ± 0.03 st at ± 0.09 syst 

3.9 dz 0.3 s tat dz 0.7 S y S t 

0.85 ± 0.02 sta t ± 0.07 syst 

Inferior conjunction (angle cut) 

2.17 ± 0.03 sta t ± 0.09 sys t 

4.1 dz 0.5 s tat dz 0.7 S yst 

1 . 1 1 ± 0.03 st at ± 0.07 syst 

Superior conjunction (angle cut) 

2. 15 ± 0.03 sta t ± 0.09 sys t 

5.0 dz 0.7 s tat dz 0.7 S y S t 

0.91 ± 0.02 stat ± 0.07 syst 

Periastron 

2. 14 ± 0.02 s tat ± 0.09 sys t 

4.1 dz 0.4 s tat dz 0.7 S y S t 

1.01 ± 0.02 sta t ± 0.07 syst 

Apastron 

2.01 ± 0.03 st at ± 0.09 sys t 

3.7 dz 0.3 s tat dz 0.7 S y S t 

0.90 ± 0.02 sta t ± 0.07 syst 



Figure 4. Overall 30 month LS I +61° 303 spectrum (red) in comparison with 
the earlier published one (black) over 8 months is shown. TeV data points taken 
by MAGIC and VERITAS are shown in gray. They are not simultaneously taken 
with the GeV data, whereas the data for the VERITAS upper limit (in green) 
are. 

(A color version of this figure is available in the online journal.) 

and corresponding best fit using the updated data set described 
in Section 2 are shown in Figure 4, together with previously 
derived results from the LAT and TeV observations. Two sets of 
TeV data are plotted: we show the non-simultaneous data points 
obtained by the Cherenkov telescope experiments MAGIC and 
VERITAS. (These data correspond to phases around 0. 6-0.7 
and represent several orbits observed in the period 2006-2008, 
before Fermi was launched.) Additionally, we show the latest 


measurements performed by VERITAS, which established a 
99% CL upper limit. 14 The new VERITAS upper limit spans 
several orbits during which, simultaneously with our LAT 
data, no detection was achieved. The LAT data along the 
whole orbit are still best described by a power law with an 
exponential cutoff. The T S value for a source emitting y-rays 
at the position of LS I +61°303 with an SED described by 
a power law with an exponential cutoff is highly significant. 
The relative T S value comparing a fit with a power law and 
a fit with a power law plus an exponential cutoff clearly 
favors the latter, at the ~20 a level. The photon index found 
is T = 2.07 ± 0.02 stat ± 0.09 S y St ; the flux above 100 MeV is 
(0.95 ± 0.01 s tat ± 0.07 syst ) x 10~ 6 photons cm~ 2 s _l , and the 
cutoff energy is 3.9 ± 0.2 stat ± 0.7 syst GeV. Results for the 
obtained T S values for each fit to different data sets are listed 
in Table 1 and all fit parameters obtained for the exponentially 
cutoff power-law models are listed in Table 2. 

Figure 4 shows that the data point at 30 GeV deviates from 
the model by more than 3er (power law with cutoff, red line). 
Although in our representation it is only one point, it is in 
itself significant, with a TS value of 67 corresponding to ~8cr. 
Therefore, and similar to the case of LS 5039, with the caveat 
of having only one point determined in the SED beyond the 
results of the fitted spectral model, we investigate the possible 
presence of a second component at HEs. As in the case for 
LS 5039, we use the likelihood ratio test to compare two 
models: model A is a power law with a cutoff and model B 
is a power law with a cutoff plus an additional power law. 
According to this test, the probability of incorrectly rejecting 
model A is 5.7 x 10“ 15 (7.8cr). The TS value for this extra 


14 We derive this differential upper limit by using the VERITAS-reported 
integral flux upper limit for phases 0.6-0. 7 (Acciari et al. 2011) assuming a 
differential spectral slope of 2.6. 
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Time [MJD] 


Figure 5. Light curve of LS I +61° 303. Each point represents one orbit of 
26.496 days, whereas the black solid line represents the three-bin-smoothed 
light curve. The black dashed line mar ks the moment of the flux change in 2009 
March. 

(A color version of this figure is available in the online journal.) 

power-law component as a whole is 172, larger than in the 
case of LS 5039, and its parameters are T = 2.5 ± 0.3 sta t and 
TL-ioGeV = (3.5 ± 0.6 st at) x 10' 10 cm -2 s' 1 . The addition of 
the HE component affects the parameters of the cutoff power 
law. The best-fit parameters for the latter, when including the 
former, are T = 2.00 ± 0.03 sta t and E cat0 ff — 2.7 ± 0.3 stat GeV 
with a flux of /L^ooMeV = (0.88 ± 0.08 stat ) x 10' 6 cm' 2 s _1 . 
Compared with the cutoff energy of 3.9 ± 0.2 GeV we obtain 
after fitting only a power law with an exponential cutoff to 
the data, the cutoff energy decreases when the additional HE 
power-law component is introduced. 

An alternative model to accommodate the deviating HE point 
in the spectrum is a fit with a broken power law. This is further 
discussed in Section 4.4. 

4.2. Light Curve 

In Figure 5, the light curve for LS I +61°303 over 30 
months is shown using orbital time bins. The black dashed 
line represents the point in time when a flux change occurred 
for LS I +61° 303; this occurred in 2009 March. In Figure 6, we 
present the power spectrum of LS I +61° 303 derived from the 
total weighted photon light curve. The power spectrum clearly 
detects the orbital flux modulation with a period of 26.71 ± 
0.05 days. This is consistent with the known orbital period of 
LS I +61°303 (Gregory 2002). The light curve clearly shows 
long-term variability. We searched for changes in the orbital 
modulation of LS I +61° 303 by dividing the aperture photometry 
light curves into six month segments and calculating the power 
spectrum of each segment, as shown in the left panel of Figure 7. 
The amplitude of the orbital modulation is estimated by fitting 
a sine wave fixed to the orbital period to each of the light curve 
segments; the results, as seen in the right panel of Figure 7, 
clearly show a decreasing trend in the orbital modulation with 
time. 

LS I +61°303 is one of the brightest sources in the y-ray 
sky and towers above all other emitters in its neighborhood. 
This allows us to compute a light curve with an orbital binning 
(26.496 days per bin) which is shown in Figure 5. Even by eye, 
it is clear that the source is highly variable on orbital timescales 
and longer. The longer term trends are evident by looking at a 
plot of the three-orbit rolling average (black line in Figure 5). 
During the first eight orbits, the flux decreases by a factor of 
~2. Then, in 2009 March, the flux appears to increase over 
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Figure 6. Lomb-Scargle periodogram of the whole LS I +61° 303 GeV data set. 
The orbital period is clearly visible. 



the course of several orbits; we take the transition point of this 
increase to be MJD 54,900. 

The flux increases significantly by 33% ± 4%, rising from a 
baseline of (0.75 ± 0.03 sta t ± 0.07 syst ) x 10~ 6 photons cm' 2 s' 1 
obtained from the first eight months of data to (1.00± 0.01 stat ± 
0.07 syst ) x 1 O' 6 photons cm' 2 s' 1 which is the average flux of 
the remaining 1.7 years of the data. Comparing the flux levels 
averaged over the same time span, eight months before and 
eight months after the flux change, we obtain a ~40% increase. 
After this flux change, the flux decreases again slowly over the 
remaining 1.7 years. The complexities of the short timescale, 
orbit-to-orbit variability make it impossible to characterize the 
exact properties of the transition from the “lower” to “higher” 
flux states. The transition likely took place over several orbits, 
however, for simplicity throughout the remainder of this analysis 
we use a transition time of MJD 54,900. 

We graphically show the flux change in Figure 8, by plotting 
the folded light curves before and after the transition in 2009 
March. The data points are folded on the Gregory (2002) period, 
with zero phase at MJD 43,366.775. Before the transition, the 
modulation was clearly seen and is compatible with the already 
published phasogram, whereas afterward, the amplitude of the 
modulation diminishes. We quantify this behavior by measuring 
the flux fraction below. Note that the data sets corresponding 
to the reported results (Abdo et al. 2009a) and what we here 
referred to as before the flux change span almost exactly the 
same time range, with the consequence of our current analysis 
essentially reproducing that previously published. The time span 
covered by our earlier publication coincidentally finished just 
prior to the onset of the flux change. The spectra derived before 
and after this flux change are shown in Figure 8, where the 
increase in flux is also obviously visible. 

4.3. Phase-resolved Spectral Analysis 

The statistics of the current data set allow us to divide the orbit 
into different phase ranges and to compute the corresponding 
spectra for different phase bins. We have divided the orbit into 
INFC and SUPC phase ranges in two different ways. First, we 
have split the LS I +61°303 orbit in two halves based on its 
geometry, as visualized, for instance, in Aragona et al. (2009). 
The SUPC phase range is defined as from phase 0.63 to phase 
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Figure 7. Left: the Lomb-Scargle periodograms of the LS I +61° 303 30 month light curve split into five consecutive segments; the earliest is at the top. The red 
dashed line indicates the known orbital period. Right: the change in the amplitude and phase of the orbital flux modulation calculated by fitting a sine wave to each of 
the light curve segments. 

(A color version of this figure is available in the online journal.) 





Energy [MeV] 


Figure 8. Left: folded light curve (100 MeV-300 GeV) of LS I +61° 303 before the flux change (blue), when the modulation is still clearly visible, compared with 
the earlier published eight month data set shown in gray. Middle: folded light curve (100 MeV-300 GeV) after the flux change, in 2009 March (red). The modulation 
gets fainter. For comparison, the previous published light curve is also plotted in gray. Right: comparison of the spectra derived before (blue) and after (red) the flux 
change in 2009 March. 

(A color version of this figure is available in the online journal.) 
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Figure 9. Spectra of LS I +61° 303 for different phase bins during the orbit are shown. Left: spectra for the geometrical cut of the orbit in superior (blue) and inferior 
(red) conjunction. Middle: spectra for the angle-based cut in superior and inferior conjunction. Right: spectra derived for periastron (red) and apastron (blue). 

(A color version of this figure is available in the online journal.) 


0.13; INFC is defined correspondingly, as the remaining half. 
We also adopted another way to separate between INFC and 
SUPC based on the angle between the compact object, the star, 
and the observer. Therefore, the orbit is not divided in two halves 
in this way, but in one piece when the compact object is in front 


of the star (corresponding to INFC): 0.244-0.507; and in another 
phase range with the same duration, centering at the exact SUPC 
phase: 0.981-0.244. As can be seen in Figure 9, the spectra 
obtained with the different cuts do not differ significantly and all 
spectral parameters are compatible within their corresponding 
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Figure 10. Modulation of the spectral parameters of LS I +61° 303 on the total 30 month light curve. For details see Section 4.3. 
(A color version of this figure is available in the online journal.) 


errors (see Table 2). We find that the flux difference between 
INFC and SUPC is of the order of 20%. For these different 
data sets, representing only portions of the orbit, we modeled 
the source with a pure power law and with a power law with 
an exponential cutoff. The ATS value comparing both fits for 
the angular based cut is 108 (INFC), which means that the 
probability of incorrectly rejecting a power law with respect to 
an exponential cutoff is 2.8 x 10 25 (10cr). For SUPC the AT S 
value is 97, which leads to a probability of 6.9 x 10~ 12 (6.8<r) 
to wrongly reject the cutoff power law. Hence, the exponential 
cutoff is preferred over the pure power law also for two parts of 
the orbit, namely, INFC or SUPC. 

We have also divided the orbit into phase ranges correspond- 
ing to periastron (half of the orbit around phase 0.275) and 
apastron (the other half, around phase 0.725), based just on the 
distance between the compact object and the star. In this case, 
neither a significant difference between the flux values for the 
two phase bins nor a difference in the spectral shape is visible. 
These results can also be seen in Figure 9. This is probably the 
result of dividing the orbit into phase ranges which contain both 
bin phases corresponding to INFC and SUPC. 

We also studied the spectral behavior of the source in phase 
bins of 0.1. For this study, we modeled LS I +61° 303 with a 
power law with a cutoff for each phase bin individually. The 
spectral parameters obtained are shown in Figure 10. Note that 
we fixed in the model the index at 2.07 to study the orbital 
behavior of the cutoff energy and we fixed the cutoff energy at 


3.9 GeV to study the behavior of the index, since both parameters 
are correlated. These values are, respectively, the results arising 
from the fit over the whole data set. A clear orbital modulation 
of the flux and spectral shape is seen. Through the periastron 
passage the spectrum gets softer and the flux is maximum, 
whereas around apastron the spectrum becomes harder and the 
flux reaches its minimum. 

4.4. Spectral Fitting 

Both the orbitally averaged spectrum and the spectra of the 
several data sets mentioned were also fitted with a broken 
power law, in addition to the pure and exponentially cutoff 
power law. All the T S values for the different fits are listed in 
Table 1. It is evident that the power law with an exponential 
cutoff or the broken power law always describe the spectral 
shape better than a pure power law does. To be precise, when 
comparing the exponentially cutoff power law with the pure 
power law, the AT S values span the range from 90 to 520; with 
the former always being statistically preferred. Fitting a broken 
power law gives almost the same results. The AT S values span 
the range from 77 to 494 and we find break energies in the 
range of 0.4-1. 7 GeV. All of the fits to the various data sets 
with an exponentially cutoff power law have slightly better T S 
values with fewer degrees of freedom than a broken power 
law, suggesting that the former model is a better description 
of the data than the latter one. The exponentially cutoff power 




The Astrophysical Journal, 749:54 (12pp), 2012 April 10 


Hadasch et al. 



0.8 


Superorbital Phase 
1.2 


c 0.6 
o 

o 

S 0.5 

U. 

x 

E 0.4 

■O 

0) 

ro 0.3 

D 

■c 

o 

s 0.2 
0.1 


- gamma-ray, Fermi 
X-ray, RXTE 


54700 54800 54900 55000 55100 55200 55300 55400 55500 55600 

MJD 


£ 

o 

.c 

Q. 

CO 

© 


1.1 


1 — 


0.8 — 


1 I 1 1 1 1 I 1 1 ' 1 I - 


average gamma-ray flux, Fermi 
average X-ray flux, RXTE 


</> 


E 

■ 1.1 “ 
o 
0 


-1 2 


- 0.8 


54700 54800 54900 55000 55100 55200 55300 55400 55500 55600 

MJD 


Figure 11. Comparison between the y - ray (blue) and the X-ray (gray) data of LS I +61° 303. Left: for each of the five separate six month periods, the 3-10 keV (gray) 
and the 100 MeV-300 GeV (blue) folded light curves are shown. Right: six month modulated flux fraction data in both energy bands, with equal color coding. 

(A color version of this figure is available in the online journal.) 


law nicely describes the curvature of the spectrum, especially 
at low energies. At HEs, above the spectral break, the broken 
power-law fits all the spectral points, even the highest one which 
might possibly be considered part of a second component. 
Statistically, we cannot distinguish which of these two fitting 
models describes best the data as a whole, but only that both of 
them are preferred over a pure power law. 

4.5. The Multi-wavelength Context 
4.5.1. X-Rays 

LS I +61° 303 has also been monitored with the 
AAT’E-Proportional Counter Array (PCA) and folded light 
curves were produced using the same ephemeris as described 
in Section 4.2. In the left panel of Figure 1 1, we show a direct 
comparison between the phasograms in X-ray and in y-rays, 
with simultaneously taken data. We divide the whole LAT data 
set into five periods of six months each and compare them with 
correspondingly obtained PCA X-ray data. The division into pe- 
riods of six months is justified in order to have enough statistics 
in y-rays for each individual time bin, such that orbit-to-orbit 
X-ray variability does not dominate the flux fraction changes. 

For RXTE- PCA, we have used the “Standard 2” mode 
for spectral analysis. Data reduction was performed using 
HEASoit 6.9. We select time intervals where the source 
elevation above Earth limb is >10° and the pointing offset 
is <0°02. PCA background light curves and spectra were 


generated using the FTOOLS task pcabackest. pcarsp was 
used to generate PCA response matrices for spectra. The 
background file used in the analysis of PCA data is the most 
recent available from the HEASARC Web site for faint sources, 
and detector breakdown events have been removed. 13 A power- 
law shape, with absorbing hydrogen column density fixed at 
0.75 x 10 22 cm" 2 (Kalberla et al. 2005; Smith et al. 2009), 
was used to fit the Standard 2 data, with the following function 
N(E) — Ke~ Na afE \E / keV) _ “, where K is a normalization at 
1 keV, a is the photoelectric cross-section, and a is the photon 
index. 

It is apparent that the X-ray modulation is always visible in 
each of these five panels, albeit with variable amplitude of flux 
modulation. We define the flux fraction as (c max — c m j n )/(c max + 
Cmin), where c max and c m i n are the maximum and minimum 
flux in the 3-10 keV found in the orbital profile analyzed 
(after background subtraction). The modulation becomes even 
stronger over time and stays stable over the last three half-year 
bins. Instead, at GeV energies, the LAT data indicate that the 
modulated fraction fades away until the variability along the 
orbit is barely visible in the last six months of our data, which is 
consistent with Figure 7. The flux fraction is plotted in the right 


15 The background file is pca_bkgd_cmfaintl7_eMv20051128.mdl and see the 
Web site: http://heasarc.gsfc.nasa.gov/docs/xte/recipes/pca_breakdown.html 
for more information on the breakdowns. The data have been bary centered 
using the FTOOLS routine faxbary using the JPL DE405 solar system 
ephemeris. 
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Figure 12. Comparison of non-contemporaneous LAT-obtained y-ray with radio data of LS I +61°303 folded with the radio super-orbital period. Top: with respect to 
the radio data compiled by Gregory (1999), mostly at 8.3 GHz. Bottom: with respect to the EW of the Ha line obtained by Zamanov et al. (1999). 

(A color version of this figure is available in the online journal.) 


top panel of Figure 1 1 . The difference between the behavior of 
the X-ray and the /-ray emission is clearly visible. In the right 
bottom panel of Figure 1 1, we show the average /- and X-ray 
fluxes fitted in each of the periods considered. We checked for 
a possible appearance of the super-orbital periodicity 1667 ± 
8 days, taking zero phase at MJD 43,366.275 and the shape of 
the outburst peak flux modulation estimated by Gregory (2002) 
in radio by considering the flux evolution and the direct count 
rate, using the full data set of the PCA observations. We find no 
evidence of the super-orbital period in X-ray or /-rays, which 
is not surprising given the short integration time in comparison 
with the super-orbital period duration. 

4.5.2. Radio and Optical 

We have also compared the LAT-obtained data with radio 
and Ha observations. For the latter, we first take into account 
the results contained in the long-term coverage presented by 
Zamanov et al. (1999). These authors have already shown that 
the equivalent width (EW) and the peak separation of the Ha 
emission line appear to vary with the super-orbital radio period 
of ~ 1600 days, this likely being the result of cyclical variations 
in the mass-loss rate of the Be companion and/or of density 
variability in the circumstellar disk. Zamanov et al. (1999) 
proposed that the variability in the EW(Ha) can be explained 
by a cyclical change in the mass-loss rate of about 25% over 
its average value. This mechanism would imply changes in the 
density of the circumstellar disk in the same range. The fact that 
the /-ray light curve would have a similar behavior to that found 
for the EW(Ha) would imply that the former is related to the 
circumstellar disk surrounding the Be star. This would naturally 
be the case if the /-rays are produced in an inter-wind shock 
formed by the collision of outflows from the compact object 
and the star itself. Dubus (2006) and Sierpowska-Bartosik & 
Torres (2009) discussed how the many uncertainties present 
in our knowledge of the LS I +61° 303 system influence the 


models for producing /-ray emission. One of the important 
parameters is the density and size (and the possible truncation) 
of the circumstellar disk. The latter influence the position of the 
inter-wind shock, and the time intervals along the orbit in which 
the compact object outflow may be balanced by the equatorial 
wind feeding it, and/or when the compact object is directly 
within the circumstellar disk itself. For recent measurements 
on very long term optical variability of Be high-mass X-ray 
binaries in the Small Magellanic Cloud, see Rajoelimanana et al. 
( 2011 ). 

We plot in Figure 12 a comparison of the LAT data, folded 
on the super-orbital radio period, with the radio data compiled 
by Gregory (1999) — most of it obtained at 8.3 GHz, using 
the Green Bank Interferometer — and the EW of the Ha line 
obtained by Zamanov et al. (1999). Care should be exercised 
when comparing our LAT data with the radio data in Zamanov 
et al. (1999), since the latter authors used the ephemeris given 
by Gregory (1999), for which the combination of the phase and 
peak flux density yielded a best-modulation period of 1584 days. 
This super-orbital period was later revised by Gregory (2002) 
to the current value of 1667 days. This latter value is the super- 
orbital period that we have used to fold the X-ray and LAT 
data, and we also use it here for the radio and Ha emission. 
The comparison is clearly non-contemporaneous, and this may 
induce problems in the interpretation of a source like LS I 
+61°303, presenting different variability phenomenology. We 
do not see any clear correlation or anti-correlation with the 
radio and Ha fluxes which may be hidden by the scarcity of 
GeV data. 

We have also compared the Fermi data with simultaneous ra- 
dio data from the Owens Valley Radio Observatory (OVRO) and 
the Arcminute Microkelvin Imager (AMI) array (Cambridge, 
UK). Since the launch of Fermi in 2008 June, the OVRO 40 m 
single-dish telescope, which is located in California (USA), has 
conducted a regular monitoring program of Galactic binaries 
(e.g., Cyg X— 3; Abdo et al. 2009c). The OVRO flux densi- 
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Figure 13. Light curve of LS I +61°303. Each point represents one orbit of 26.496 days. The black dashed line marks the moment of the flux change in 2009 March, 
together with the simultaneously taken radio data by the AMI and OVRO instruments (both operating at 15 GHz) are shown in gray. 

(A color version of this figure is available in the online journal.) 


ties are measured in a single 3 GHz wide band centered on 
15 GHz. A complete description of the OVRO 40 m telescope 
and calibration strategy can be found in Richards et al. (201 1). 
Furthermore, we also use complementary observations provided 
by AMI, consisting of a set of eight 13m antennas with a max- 
imum baseline of ~120 m. AMI observations are conducted 
with a 6 GHz bandwidth receivers also centered at 15 GHz. See 
Zwart et al. (2008) for more details on the AMI interferometer, 
which is mostly used for the study of the cosmic microwave 
background. By folding these radio data, no super-orbital mod- 
ulation could be seen, which could be due to the poor coverage 
of a whole super-orbit. A direct comparison of the GeV and the 
radio data is shown in Figure 13 where the flux over time is 
plotted. No correlation of the data points can be seen either. In 
summary, we did not find any correlation between the GeV and 
the radio band, either in archived or in simultaneous data. 

5. CONCLUDING REMARKS 

After analyzing a data set comprising 2.5 years of Fermi-LXT 
observations of the two binaries LS 5039 and LS I +61°303, we 
note several changes with respect to the initial reports of Abdo 
et al. (2009a, 2009b). These were produced either because the 
accumulation of a longer observation time allowed us to make 
distinctions that were earlier impossible (valid for both sources), 
or because the behavior of the source changed (valid for 
LS I +61°303). On one hand, the statistics are now sufficient 
to divide the data set of both sources in INFC and SUPC and 
to show that a power law with a cutoff describes the spectra 
obtained in both conjunctions better than a pure power law. The 
cutoff is similar to that found in the many other GeV pulsars 
discovered by hermi-LKY. However, we have found that both 
LS 5039 and LS I +61° 303 show an excess of the HE GeV 
emission beyond what is expected from an exponentially cutoff 
power law. While the HE data are significantly in excess of the 
exponentially cutoff power law, there are insufficient statistics 
at these HEs to model the excess with an additional spectral 
component. 


The process(es) generating such a component in the cases of 
LS 5039 and LS I +61°303 is unclear and may even be different 
in the two sources. However, such a second component would 
present a possible connection between the GeV and TeV spectra 
in both sources. Collecting more data, and therefore more 
statistics, will allow us to prove or discard it in the future. The 
lack of data points at HEs also affects, particularly in the case of 
LS I +61° 303, the distinction between an exponentially cut and 
a broken power law. Currently, both are certainly preferred over 
a pure power law, but differences in the significances provided 
between the former are minor. 

We have noted that, whereas LS 5039 shows stable emission 
over time and also a stable orbital modulation, LS I +61° 303 
shows a change in flux in 2009 March. Afterward, the or- 
bital modulation decreases (see the bottommost left panels of 
Figure 11) and the orbital period could not be detected in the 
GeV data. LS I +61° 303 has also presented a complex, concur- 
rent behavior at higher energies. At TeV, for approximately the 
last two years, it seems to have been in a low state in comparison 
with the flux level that led to its discovery (Aleksic et al. 2011; 
Acciari et al. 2011). Additionally, it was detected once — after 
4.2 hr of observations — near the periastron, where the system 
was never seen at TeV energies before (Acciari et al. 2011). 
Both of these aspects of the LS I +61° 303 phenomenology — as 
well as the phase location of the TeV maxima in general — would 
require modifications of simple inverse Compton models in the 
scenarios usually put forward for the source. The idea of a mag- 
netar compact object in LS I +61°303 is suggested by Torres 
et al. (2012) to explain the changing TeV behavior of this object 
and could potentially describe the diminishing orbital modu- 
lation now observed at GeV fluxes. However, as explained by 
Torres et al. (2012), detailed numerical simulations are required 
to verify if this scenario can accurately reproduce the observed 
GeV and TeV emission. 

The lower-energy, multi-wavelength picture of LS I +61° 303 
is still unclear, with the orbital modulation in GeV (X-ray) fading 
(increasing) with time. Because the GeV data do not cover 
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a whole super-orbital period, a conclusion about the possible 
relationship of the GeV emission with the super-orbital behavior 
of LS I +61° 303 cannot be drawn yet. Continued monitoring of 
the source with the Fermi - LAT will allow for the study of such 
long-term cycles, if any. Zamanov et al. (1999) concluded that 
if the ~4 year modulation in radio and Ha is due to changes in 
the circumstellar disk density, then it could be detected in X- or 
y-rays as well, since the HE emission of a putative neutron star 
depends on the density of the surrounding matter with which 
it interacts. Indeed, this modulation was recently revealed in 
X-rays (Li et al. 2012). 

The nature of these systems is thus still unclear and multi- 
wavelength observations, including Fermi-LXT monitoring, as 
well as further theoretical studies are expected to bring new 
insights. 
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